#ifndef CHANNELFLOW_HPP
#define CHANNELFLOW_HPP

#include "edgeTraversal.hpp"
#include "sweSolver.hpp"
#include "solverBase.hpp"
#include "meshTypes.hpp"
#include "meshTypeFunctions.hpp"
#include "magic_enum.hpp"
#include "vtkMeshExport.hpp"
#include "bc_freeOutflow_velocity.hpp"
#include "bc_freeOutflow_waterdepth.hpp"
#include "bc_constantWaterDepth_velocity.hpp"
#include "bc_staticValue.hpp"
#include "bc_zeroGradient.hpp"
#include "bc_slip_velocity.hpp"

using namespace Eigen;
using namespace hms;

struct channelFlow{
	template<typename MeshType>
	void runCase(){
		// mesh properties
		Index 
			nx {10},
			ny {10},
			ng {1};
		scalar dx {0.5};
		// mesh generation
		auto mesh { uniMeshAs<MeshType>(nx,ny,dx,dx) };
		addStandardBoundaryPatches(mesh);

		bool useFriction {true};


		// solver properties
		SWESolver<MeshType> solver;
		solver.edgeTraversal(EdgeTraversal::blockwise_secondOrderAlt);
		solver.useFriction(useFriction);
		solver.useSourceFields(false);
		solver.useMuscl(true);
		solver.limiterFunction(LimiterFunction::minMod);
		solver.maxCourant(0.3);

		// initial values and boundary conditions
		solver.v() = Field{mesh, 2}; // velocity
		solver.d() = Field{mesh, 1}; // depth
		solver.z() = Field{mesh, 1}; // bed elevation
		solver.v().values().colwise() = Array2s{0.2, 0.0};
		solver.z().values() = 0.0;

		solver.d().values() = 1.0; // left half of the domain is d=1.5, right half d=1.0

		solver.d().values().reshaped((nx+2*ng), (ny+2*ng)).transpose().leftCols(0.5*(nx+2*ng)) = MatrixXd::Constant(ny+2, 0.5*(nx+2), 1.5);

		solver.n() = ConstantField<Uniform, MeshType>{ mesh, ArrayNs<1>{0.025} };

		const BoundaryPatch& left   { mesh.boundaryPatch("left") };
		const BoundaryPatch& right  { mesh.boundaryPatch("right") };
		const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
		const BoundaryPatch& top    { mesh.boundaryPatch("top") };

		solver.d().setBoundaryCondition( left, StaticValue{}, Array<scalar, 1, 1>{1.5}, true);
		solver.d().setBoundaryCondition( right, FreeOutflow_waterdepth{});
		solver.d().setBoundaryCondition( top, ZeroGradient{} );
		solver.d().setBoundaryCondition( bottom, ZeroGradient{} );

		solver.v().setBoundaryCondition( left, ConstantWaterDepth_velocity{}, solver.d() );
		solver.v().setBoundaryCondition( right, FreeOutflow_velocity{} );
		solver.v().setBoundaryCondition( top, Slip_velocity{} );
		solver.v().setBoundaryCondition( bottom, Slip_velocity{} );

		solver.z().setBoundaryCondition( left, ZeroGradient{} );
		solver.z().setBoundaryCondition( right, ZeroGradient{} );
		solver.z().setBoundaryCondition( top, ZeroGradient{} );
		solver.z().setBoundaryCondition( bottom, ZeroGradient{} );


		solver.initialise();

		// simulation time
		Time time {120}; // 10 seconds of real time.
		std::cout << std::setprecision(5) << std::fixed; // digits afte comma = 5
		std::cout << "Time: " << time.current() << "s";

		// Timer for benchmarking purposes
		Timer t;

		// go simulate 10 seconds of real time
		while( !time.atEnd() ){
			solver.advance( time );
			std::cout << "\b\b\b\b\b\b\b\b";
			std::cout << time.current() << "s";
		}

		// print information about the performance time.
		double elapsed { t.elapsed() }; 
		double execTimeSum { 0.0 };

		for ( auto execTime : solver.execTimings() )
			execTimeSum += execTime;

		std::cout
			<< "\nComputation time: " << elapsed << "s"
			<< "\nSum of execution times logged by solver: " << execTimeSum << "s"
			<< "\n\t" << 100 * solver.execTimings()[0] / execTimeSum
				<< "% for time step calculation"
				<< " (" << solver.execTimings()[0] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[1] / execTimeSum
				<< "% for fluxes (including slope)"
				<< " (" << solver.execTimings()[1] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[2] / execTimeSum
				<< "% for sources " << ( useFriction ?
					std::string{ "(including friction)" } :
					std::string{ "(friction disabled)" } )
				<< " (" << solver.execTimings()[2] << "s)"	
			<< "\n\t" << 100 * solver.execTimings()[3] / execTimeSum
				<< "% for new state"
				<< " (" << solver.execTimings()[3] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[4] / execTimeSum
				<< "% for updating boundaries"
				<< " (" << solver.execTimings()[4] << "s)"
			<< "\n";
	}
};

#endif